Brain organization into resting state networks emerges from the connectome at 

criticality 
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The relation between large-scale brain structure and function is an outstanding open problem 
in neuroscience. We approach this problem by studying the dynamical regime under which re- 
alistic spatio-temporal patterns of brain activity emerge from the empirically derived network of 
human brain neuroanatomical connections. The results show that critical dynamics unfolding on 
the structural connectivity of the human brain allow the recovery of many key experimental findings 
obtained with functional Magnetic Resonance Imaging (fMRI), such as divergence of the correlation 
length, anomalous scaling of correlation fluctuations, and the emergence of large-scale resting state 
networks. 



Understanding the relation between brain architecture 
and the behavioral repertoire it produces is a central 
question in neuroscience. In that direction, important 
efforts over recent years have been devoted to map the 
large-scale structure of the human cortex, including at- 
tempts to build brain structural connectivity matrices 
from imaging data. An example is the connectivity ma- 
trix of the entire human brain, recently derived from fiber 
densities measured between a large number (500-4000) 
of homogeneously distributed brain regions [1]. This and 
related work encompasses a large collaborative project 
dubbed the brain "connectome" [2], whose ultimate goal 
is to understand in detail the architecture of whole-brain 
connectivity. The question discussed in this Letter is 
how much of the experimentally observed brain dynam- 
ics can be then expounded from such detailed knowledge 
of brain connectivity. In other words, and using a simple 
analogy: how much of the large-scale traffic dynamics of 
a city can be predicted from a very detailed map of all 
streets?. The results presented in this Letter show that 
very relevant aspects of brain dynamics can be predicted 
from the structure providing that the underlying dynam- 
ics are critical. 

To guide our comparison with available experimental 
results, we choose to concentrate on robust findings con- 
cerning brain dynamics. Specifically, we wish to ask how 
spontaneous brain dynamics at the large scale organize 
into the relatively few spatio-temporal patterns revealed 
experimentally in recent years [3j. This is important be- 
cause a wide range of experiments using functional Mag- 
netic Resonance Imaging (fMRI) have repeatedly em- 
phasized that these spatial clusters of coherent activity, 
termed Resting State Networks (RSN) [4| , are associated 
with specific sensory, cognitive and behavioral functions 



[5[. Their functional role is currently being investigated 
both in health and disease under a variety of angles. 
Of interest here are studies showing that the RSN activity 
exhibits peculiar scaling properties, resembling dynamics 
near the critical point of a second order phase transition 
@-@|. In addition, there is evidence showing that the 
brain at rest spends most of the time wandering near 
a critical point [9[. These empirical findings are consis- 
tent with results obtained with computational modeling 



Here we study whether a simple dynamical model run- 
ning over the empirical structure of neuroanatomical con- 
nections [1] suffices to replicate the aforementioned fun- 
damental features of spontaneous brain activity repeat- 
edly seen in fMRI experiments. The model consists of a 
network of interconnected nodes (i.e., the connectome) , 
together with a dynamical rule. The matrix of connec- 
tions follows the neuroanatomical connectivity described 
recently by Hagmann et al. [1], who studied healthy hu- 
man subjects and reported the average fiber tract den- 
sity between any two brain areas (from a gray matter 
parcellation into 998 areas). To complete the model we 
need to specify the dynamics of each node. For simplic- 
ity, we choose discrete state excitable dynamics following 
the Greenberg-Hastings model [13]. Thus, each node is 
assigned a state which can be one of three -quiescent Q, 
excited E, or refractory R - with the dynamics deter- 
mined by the transition rules: 1) Q —> E with a small 
probability n (~ 10 -3 ), or if the sum of the connection 
weights Wij with the active neighbors (j) is higher than 
a threshold T, i.e., ^2 Wij > T and Q — >• Q otherwise; 
2) E — >• R always; 3) R — >• Q with a small probability 
V2 (~ 10 _1 ) delaying the transition from the R to the Q 
state for some time steps. We held fixed parameters r\ 
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FIG. 1. The model. Panel A shows the model adjacency 
matrix (i.e., the connectome) and Panel B the edge weights 
distribution. Data obtained from the white matter tracts con- 
necting a parcellation of 998 nodes covering the entire hu- 
man brain [l|]. RH and LH refers to the right and left brain 
hemisphere, respectively. Panel C shows the giant cluster's 
size, (i.e., the order parameter, Si, solid line) and the second 
largest cluster's size (&, dashed line) as a function of the 
threshold T (control parameter) as well as the critical point 
T c ~ 0.05. Panel D illustrates examples of clusters (denoted 
with different colors) at the three T values denoted with col- 
ored markers in Panel C. 



and T2 , which determine the time scales of self-excitation 
and of recovery from the excited state, respectively, and 
changed T. For the numerical analyses, the time series 
of each node was binarized by assigning state E = 1 and 
the remaining states into 0's and convolved with a stan- 
dard hemodynamic response function [14] which mimics 
the coupling between neural and metabolic activity mea- 
sured in the fMRI experiments. 

As depicted in Figure 1, the dynamics of the model 
show a transition as a function of the threshold T. For 
relatively small values of T even the weakest connec- 
tions are enough for the activity to spread, resulting in 
a regime with a relatively high activity level. On the 
contrary, for high values of T the activity only flows 
through the few strongest connections and therefore the 
overall activity decreases. To characterize the transition 
between these regimes an order parameter was defined 
considering the sizes of the active clusters. Clusters are 
groups of nodes simultaneously activated and linked to 
each other through a non zero w. At each time step the 



size of the largest cluster (Si) and the second largest clus- 
ter (S2) were computed. These calculations (see Panel C 
of Fig. [I]) unveil a transition between a phase in which a 
giant cluster covers ~ 15% of the system (while the sec- 
ond largest cluster is of negligible size) and another phase 
in which only scarce activations occur and the nodes fail 
to coalesce into large clusters. At an intermediate value, 
a critical point (T c in Fig. 1C) can be identified by the 
peak in the size of the second largest cluster, as done 
usually in percolation [15] as well as recently in human 
fMRI experiments Q. 

We compare now the dynamics of the model with previ- 
ous experimental results. In particular, we are interested 
in two robust features exhibited by the spontaneous ac- 
tivity of human brain RSN [16]: 1) the correlation length 
of brain activity diverges with size, and 2) the variance 
of the short-term correlations between pairs of brain sites 
remains high, independently of the number of pairs con- 
sidered. Since these two properties are often seen as 
generic features of criticality, we decided to explore first 
whether the model exhibits similar dynamics. To be able 
to compare with the experimental results, we identified 
each node of the model as belonging to the closest human 
RSN. This was done by matching the node's coordinates 
provided by Hagmann et al [lfwith a spatial mask of the 
RSN [4] as done previously [16|. Nodes that were far- 
ther than 1 cm away from the closest RSN cluster were 
discarded from the analysis. 
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FIG. 2. The correlation length £ of the activity in the model 
near T c diverges with the cluster size, as reported for human 
brain data [lq |. Panel A shows the correlation function C(r) 
computed from human data (Exp.) and from the model at 
T c (colored lines are used for the different clusters). The cor- 
relation length £ is the distance r where C(r) = 0, seen here 
to span a range (denoted with the arrows). Panel B shows 
the £ values for the functions plotted on panel A, demonstrat- 
ing that £ ~ N 1 / 3 (dashed line), both in the experiment and 
model data. 



Divergence of the correlation length. The correla- 
tion length represents the average distance at which two 
points in the system behave independently, and is known 
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to diverge at criticality [17j . Following a standard proce- 
dure 16, [ItJ , we computed the average correlation func- 
tion of the signal fluctuations between all pairs of nodes in 
each cluster which are separated by a distance r, yielding 
the correlation function C(r) (details on this computa- 
tion can be found in the Supp. Info, as well as in |16|). 

Fig. 2A corresponds to the two point correlation func- 
tion as a function of distance for all cluster sizes obtained 
experimentally and in the model at T c . Panel B shows 
the dependence of £ with the cluster size N. Both numer- 
ical estimations clearly show that near T c the diver gen ce 
of the correlation length found in the experiments [16| is 
reproduced by the model. 
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FIG. 3. The short-term correlation < C > in the model at 
T c exhibits transient fluctuations at all cluster sizes as seen in 
human brain data [lfj]. Panel A: Examples of < C > fluctua- 
tions for three cluster sizes N at critical (T c ) and subcritical 
(T < T c ) dynamical regimes of the model, as well as for the 
human brain data (Exp.). Panel B shows that the variance 
of the fluctuations in < C > remains approximately constant 
only for the human brain data (empty black circles) and for 
the model at T c (filled red circles). For T < T c (filled blue 
circles) the variance decreases faster with N. The three small 
arrows denote the sizes used in the examples of Panel A. Panel 
C shows a plot of the distance between the scaling of the fluc- 
tuations of the human fMRI data and those from the model 
for a wide range of T. Notice that the best agreement occurs 
for T c . Colored markers in Panel C correspond to the T c and 
T < T c . values used in Panels A and B. 

Temporal fluctuations of the mean correlation in the 
RSN. As recently shown [16|, the time evolution of the 
correlation within these patterns exhibits bursts of high 
correlation intermixed with instances of dis-coordination. 
Panel A in Fig. 3 shows examples of the fluctuations in 
the short-term mean correlation < C > between all pairs 
of nodes within a given cluster, in this case calculated in 
non-overlapping time windows of 10 steps. At the critical 



state, the variance of < C > is of the same order as ob- 
served in the experiments for different cluster sizes. Panel 
B shows the dependence of the fluctuations in < C > 
with the cluster size N. At the subcritical regime the 
fluctuations decrease as , which reveals the asynchrony 
of the active nodes. On the other hand, at the critical 
state they remain constant, similar to what is observed in 
experimental fMRI data [16| . In the supercritical regime 
the vanishinly low level of activity prevents high correla- 
tions and therefore the fluctuations amplitude are close 
to zero for all cluster sizes. To compare with the human 
fMRI experimental data, we computed the root mean 
square distance between the model (m) and the experi- 



mental data (e) as A 
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where the 



sum is over all clusters N c . Panel C shows that the dis- 
tance is minimum precisely at the critical point of the 
model. 

RSN spatial patterns emerge only at criticality. One 
of the most revealing features of large-scale spontaneous 
brain activity is its spatial organization into RSN. Their 
functional relevance is highlighted by the fact that the 
spontaneous activity closely parallels brain activation 
patterns seen during task execution [5]. We studied 
whether these patterns can be seen in the model using 
the same methods employed to reveal RSN in experi- 
mental data (Independent Component Analysis -ICA) as 
implemented in the FSL MELODIC software 0. Before 
proceeding to that, we needed to identify the cortical lo- 
cations of the model's 998 nodes. For that purpose, the 
brain gray matter was parcellated via a random growth 
algorithm (see Supp. Info) using the coordinates of the 
998 nodes as seeds. After that, the model time series of 
each region of interest was assigned to the corresponding 
parcellation patch (plus Gaussian noise of 0.15 times the 
variance of the signal), then an ICA decomposition into 
8 independent components was performed (a total of 100 
trials were performed for each T value). For each compo- 
nent, we computed the maximal correlation of the model 
with the spatial location of a set of well-established hu- 
man RSN [4] obtained with fMRI. As seen in Fig. 4, the 
correlation < C e - m > peaks near T c in all cases. This 
high correlation implies that the RSN coordinated spon- 
taneous activity unfolds at the same anatomical locations 
both in the human brain and in the model close to T c , 
something already evident by visual inspection of many 
of the patterns presented in Fig. 4. 

Discussion. To the best of our knowledge, this is the 
first demonstration that a hybrid modeling approach (re- 
alistic anatomical connectivity plus a simple dynamical 
rule) suffices to capture relevant spatio-temporal aspects 
of dynamics, provided that the dynamical regime is crit- 
ical. These aspects include generic features of critical 
systems, but also the emergence of structures having 
a well-established neurobiological meaning, namely the 
cortical RSN. While experimental evidence for the afore- 
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FIG. 4. The spatial organization of the human brain RSN 
emerges spontaneously in the model near T c . RSN obtained 
from the model with ICA are shown for three values of T: 
T < T c (0.03) , T > T c (0.1) and the correlation maxima 
with the human RSN from [4] ("fMRI"), which is always ~ 
T c (red marker). Right column shows the mean correlation 
< Ce-m > between the spatial location of the experimental 
human RSN and those obtained from the model as a function 
of T. Abbreviations: Medial visual (VisM), lateral visual 
(VisL), auditory (Aud), sensory-motor (SM), default mode 
network (DMN), executive control (EC), dorsal visual stream 
left (DorL) and right (DorR). Numbers beneath each brain 
slice denote its horizontal coordinate. 



mentioned signatures of criticality in brain systems was 
already discussed [7], in this Letter we made a stronger 
point: in the model, the critical regime appears as a nec- 
essary condition for the emergence of neurobiologically 
relevant aspects of brain dynamics. Our result also rep- 
resents an important first step in the direction of realistic 
hybrid computational modeling of large-scale brain func- 
tion both in health and disease. As an example, many al- 
tered brain states are associated with RSN alterations, a 



prominent example being the loss of consciousness in the 
comatose state [19[. In light of the present results such 
brain state alterations can be regarded as a displacement 
from an optimal dynamical point. 

Summarizing, the results show that by endowing with 
critical dynamics the brain network of anatomical con- 
nections (or connectome), key observations about brain 
dynamics can be replicated. These results contribute to 
close the gap between structural and functional network 
connectivity in the human brain, by emphasizing the dy- 
namical regime at which models should predict a wide 
range of observations about large scale brain function. 
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